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We theoretically show that an interlayer bias voltage in the AB-stacked bilayer graphene nanorib¬ 
bons with armchair edges induces an electric polarization along the ribbon. Both tight-binding 
and ab initio calculations consistently indicate that when the bias voltage is weak, the polarization 
shows opposite signs depending on the ribbon width modulo three. This nontrivial dependence is 
explained using a two-band effective model. A strong limit of the bias voltage in the tight-binding 
model shows either one-third or zero polarization, which agrees with topological argument. 
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Monolayer graphene nanoribbons (GNRs) show vari¬ 
ous energy bands depending on the edge orientation and 
the width of the nanoribbons [i|,Q. When the GNRs have 
armchair edges, the energy bands become gapped or gap¬ 
less, depending on the width. Like monolayer GNRs, the 
energy bands in the AB-stacked bilayer GNRs (Fig. [1]) 
also depend sensitively on the edges. Namely, the zigzag 
bilayer GNRs show localized edge states Q , whereas the 
armchair bilayer GNRs vary from insulator {N = 31 or 
31 -\- 1, 1: integer) to metal {N = 31 -\- 2) by chang¬ 
ing the width N in a. tight-binding (TB) model (Fig. [T]) 
[ij. In addition, external fields play scientifically and 
technologically important roles in atomic-layer materi¬ 
als, such as bilayer graphene. The external electric field 
opens up the fundamental gap in the AB-stacked bilayer 
graphene ii which otherwise possesses massive and 
gapless parabolic bands in the low energy region Q. 

In this Letter, we theoretically show that an external 
interlayer bias voltage in the AB-stacked bilayer GNR 
with armchair edges induces a polarization along the rib¬ 
bon direction. We use two methods: calculation on the 
TB model and ab initio calculations. Both two methods 
consistently show that when the bias voltage is weak, 
the polarization shows a nontrivial dependence on the 
ribbon width, having opposite signs depending on the 
width modulo three. A strong limit of the bias voltage 
shows either one-third or zero polarization in the unit of 
the electron charge, which agrees with topological argu¬ 
ment. We then discuss that the present theory applies to 
a wide variety of atomic-layer compounds. Thus nanos¬ 
tructure which breaks bulk symmetries allows novel re¬ 
sponses which are absent in the bulk. 

We first discuss symmetry requirement for the trans¬ 
verse response of the polarization along the ribbon in¬ 
duced by the interlayer bias. We take the x and z axis 
in the direction normal to the edges and the bilayer, re¬ 
spectively, and the y axis along the ribbon (Fig. [T]). For 
armchair or chiral edges, when the interlayer bias volt¬ 
age is zero, inversion and C 2 x symmetries are preserved, 



FIG. 1. (Color online) Structure of the AB-stacked bilayer 
GNR with the armchair edges. (a),(b) The black (blue) lines 
represents the bondings in the lower (upper) layer, (a) N is 
the width corresponding to a number of rows. The unit cell of 
the ribbons (red dashed-line box) contains AN atoms. nA(A’) 
and nB(B’) represent the sublattice in the lower (upper) layer 
within nth row. (b) The green arrows denote the interlayer 
hoppings between the nB and nA’ sites forming “dimers”. 

which prohibit emergence of polarization. Interlayer volt¬ 
age breaks both symmetries, resulting in a polarization 
along the ribbon. For zigzag edges, x^-plane mirror sym¬ 
metry is preserved in addition to inversion symmetry, and 
it prohibits emergence of polarization along the ribbon 
(y) direction. Because the interlayer voltage does not 
break this mirror symmetry, it does not induce polariza¬ 
tion for zigzag ribbons. 

First we numerically calculate the polarization for a 
spinless TB model. 

H = X! (1) 

The first term describes the hoppings with the amplitude 
tij, for which we only consider the nearest-neighbor in¬ 
tralayer hopping t and the interlayer hopping t± within 
a “dimer”. Here, we set t± = 0.13t {t > 0) according to 
Ref[3- The second term represents the interlayer bias U 
and takes -|-1 (—1) for the upper (lower) layers. 
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From the TB model, we calculate the electronic con¬ 
tribution of the polarization P in terms of the Berry con¬ 
nection within the modern theory of polarization [H. 
It is calculated as a change of polarization P{U) — P{0) 
by changing the interlayer bias voltage U. This calcula¬ 
tion works only for insulators, and therefore we restrict 
ourselves to the insulating GNRs whose width is N = 31 
or 3/ -I- 1 {I- integer). Because P(0) = 0 by inversion 
symmetry, we obtain P{U) numerically. 

Using this method, we find that in the bilayer GNRs 
with the armchair edges, the polarization arises for 
nonzero interlayer bias voltage U. Figures [2])a)(b) are 
our numerical results for various widths N. In reality, 
feasible values of U may be limited to about |f/| < 0.15t; 
nevertheless we show the results for much larger U in the 
figure, to show consistency for the large U limit. Inter¬ 
estingly, the behavior of the polarization is classified into 
two classes, N = 31 and N = 31-1-1. For N = 31, the po¬ 
larization goes to zero for U ±oo, while for iV = 31 -I-1 
the polarization goes to ±e/3, where —e (e > 0) repre¬ 
sents the electron charge. Furthermore, the slope around 
[/ ~ 0 has opposite signs between the two classes. The 
slope at [/ ^ 0 is steeper for wider ribbons. In the in¬ 
termediate range of U the polarization oscillates as U is 
changed. This oscillation accompanies a change of band 
structure around the Fermi energy, formed by a number 
of minibands from a finite-size effect. 



u/t u/t 

FIG. 2. (Color online) Numerical results of the polarization 
from the TB model in response to the interlayer bias. The 
interlayer hopping parameter t± — 0.13t is fixed, (a) N = 31 
{N = 27,39) and (b) N = 31 + 1 (N = 28,40). The results 
for N = 39, 40 (green) are vertically offset by 0.2. 

The dependence of the asymptotic behavior at 17 —>• 
±oo on the ribbon width can be physically understood 
as follows. Because there are two electrons per row in 
the unit cell, two electrons lie on the lower layer under 
the strong interlayer bias voltage U ^ -boo. As a result, 
compared from U ^ 0, the two electrons are displaced 
by a/6 on average, and therefore the polarization is P = 
(a/6) • 2Ne/a = Ne/3 (mod e) per unit length (Fig. [3|), 
where a is the lattice constant. Here the polarization is 
defined modulo e Hence we obtain P = 0 (mod e) 

and P = e/3 (mod e) for N = 31 and 31-1-1, respectively. 
It totally agrees with our numerical calculations. 

Next, we focus on the region of the weak interlayer bias 
voltage in Fig. [2] To understand the novel behavior of 



FIG. 3. (Color online) Emergence of the electric dipole mo¬ 
ment in the strong limit of the interlayer bias voltage U. (a) 
and (b) show top and side views of a part of the nanoribbon. 
0 and —2e represent the charge at each layer for large U. The 
dashed and solid red arrows represent the whole and the y 
component of the dipole moment. 


the slope, we construct a simple two-band (2B) effective 
Hamiltonian for the weakly biased GNRs, by retaining 
only the highest occupied band and the lowest unoccu¬ 
pied band. To this end, we begin with the analytic forms 
of the eigenstates of the TB model at fc = 0 and U = 0 
[l|. The eigenvalue equation at fc = 0, C/ = 0 is written 


( ° 

\^t(2cos0 + 1) 


t{2cos6 -b 1) 
±t± 



where corresponds to the sum and the difference be¬ 
tween the amplitudes at the A and B’ sublattices, respec¬ 
tively, and is defined similarly for the B and A’ sub¬ 
lattices (see the Supplemental Material [Uj for details). 
0 is a phase difference of an electronic wave between the 
neighboring rows, forming a standing wave in the ribbon. 
Its eigenvalues are 


= ± Y + 7^ (y) +t2(2cos6i-b 1)2, (3) 

where q = ztl. From the boundary condition, we get 
9 = 9^ = ]v^7r, r = 1, 2,..., iV. Here, response of the 
polarization P to an external perturbation is given by 
the Berry curvature [^. 13. [l^ . Therefore, the eigenstates 
close to fc = 0, where the band structure has a direct gap 
when U 0, contributes considerably to the polarization. 
Hence, from the analytic forms of the eigenstates of the 
TB model at fc = 0 and P = 0 [ij , we retain only the low¬ 
est unoccupied state |-b) and the highest occupied state 
|—). Their energy eigenvalues are given by ±(? 0 i where 

go = —t±/2+d, d = ^{t±/2Y -b P(2 cosd^^^^ -b 1)^, and 

^ 21+1 ~ 7r(2Z-bl)/(iV-bl) for both N = 31 and A^ = 3/-bl. 
By using these two eigenstates, we construct a 2B model 
which describes the energy bands around the Fermi en¬ 
ergy when fc ^ 0 and U The 2B Hamiltonian to 

the first order in fc and U is derived as 


Heff = hiU ax + h2kay + goa^, (4) 


where hi = t±/{Ad), h 2 = af^(cos— 1)(2 cosd^.,.^ -b 
l)/(3d) and a^^y^z are the Pauli matrices for the space 
spanned by the eigenstates |±) [Ilj. We note that width 
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FIG. 4. (Color online) Comparison between the results from 
the TB model (red solid line) and the 2B model (green dashed 
line), (a) N = 31 {N = 27,39) and (b) = 3/ + 1 (IV = 

28,40). The results for N = 39, 40 are vertically offset by 0.1. 


dependence appears through From the 2B Hamil¬ 

tonian, we calculate the polarization P{U) for small U, 
as shown in Fig|4l They well agree with the results of the 
TB model in the U ^ 0 region, including the sign of the 
slope, for N > 20. For N < 20, the gap in the TB model 
is non-monotonous as a function of the interlayer bias U ; 
for small U the gap decreases as a function oi U. It is 
not reproduced in the 2B model where the gap always 
increases with U. This leads to differences between the 
two models. 

The width dependence is understood from the analytic 
formula for the slope of P(U)\ 


dP ie 
dU 2 tt 


dk- 


{Ek- - Ek+Y 


Here Ekj and \ukj) are the eigenvalue and the eigenstate 
of the 2B Hamiltonian for the j band. Therefore, the 
slope at [/ = 0 is 


— e ^ I, 1 

dU u=o 2a ^ ^ go\/90 + (/i27r/a)2’ 


( 6 ) 


and its sign is given by —sgn(/i 2 ). From 2tt/3 — = 

=F7r/(37V -I- 3) for N = 31 and 3Z -I- 1, the sign of the slope 
is negative for N = 31 and positive for N = 31 + 1, in 
agreement with the results for various widths of GNRs. 
Asymptotic form for a wider ribbon is evaluated as 


dP 


u=o 


3tj e 


{N + lf. 


(7) 


where the signs ± is — for = 31 and -|- for iV = 3Z -I- 1. 


These behaviors are confirmed by ab initio calcula¬ 
tions based on density functional theory (DFT). We 
perform the electronic structure calculation of hydrogen 
terminated AB-stacked bilayer GNRs within the local- 
density approximation (LDA) [l^, [14 [based on DFT us¬ 
ing Quantum Espresso package 1^. We use ultra- 
soft pseudopotentials [l^ and plane-wave basis sets to de¬ 
scribe the charge densities and wave functions with cutoff 
energies of 30Ry and 300Ry, respectively. The supercell 
approach is used and the distances of neighboring bilayer 


GNRs along the i-axis and the z-axis are at least 10 and 
30 A, respectively. The geometries are fully optimized. 
To discuss the effect of the external electric field, we ap¬ 
ply a periodic zigzag potential along the z-axis in the 
supercell. Under the external field, E, we obtain the 
band structure with 48 x 1 x 1 k-points and calculate the 
electric polarization in terms of the Berry connection. 

The dot symbols in Fig. [5] (a),(b) represent the DFT re¬ 
sults on polarization for hydrogen terminated AB-stacked 
bilayer GNRs under the electric field E. Here, we put 
P{E = 0) = 0 by symmetry. Apparently, the signs 
of the slopes of P obtained by DFT calculations com¬ 
pletely agree with those from the TB and the 2B model. 
Furthermore, the N dependence of the polarization, i.e., 
dP{U)/dU\^^g oc (A^-l-1)^ in Eq. JT]), is well reproduced 
in the DFT results (see Fig. [5](c)). These indicate that 
the simple TB model, and consequently its 2B effective 
model, well capture the key features of the polarization in 
this system. To compare the polarization values with the 
TB model quantitatively, we relate the electric field, E to 
the on-site energies U in the TB model. To this end, we 
construct maximally localized Wannier functions for the 
valence bands using carbon tr and tt orbitals 17 , and 

the result is shown in Fig. [5] (d). The obtained on-site 
energies for tt orbitals in each layer are almost indepen¬ 
dent of the position of the orbitals except at the edges. 
Thus, in Fig. [5] (d), we use average values excluding the 
edges. 

As shown in Fig. [5] (d), we find that for weak elec¬ 
tric field \E\ < 0.2V/A, E and U are almost linear, 
whereas their proportionality constant depends on the 
ribbon width. Using this correspondence the results on 
the TB model over various U is translated into the de¬ 
pendence on the electric field E, as shown as the dotted 
lines in Fig.[S](a) and (b). We notice that the result from 
DFT and that from the TB models have similar tenden¬ 
cies, whereas they are different by a factor of two smaller 
or larger, depending on the series N = 31 and N = 31 + 1. 
This difference between DFT and the TB model can be 
partly attributed to the difference of the gap size. In 
the result of the 2B model in Eq. ©, the polarization is 
inversely proportional to the gap size, because h27r/a is 
much larger than go for the given parameters. Actually, 
the gap size obtained from DFT is smaller (larger) than 
that of the TB model for N = 31 {N = 3/ -I- 1) even at 
U = 0. Therefore, in order to incorporate this difference 
of the gap size, we rescale the results of the polarization 
of the TB model by the ratio between the gaps from the 
DFT and that of the TB model. This rescaling enhances 
(suppresses) the polarization ioi N = 31 {N = 31 + 1). 
After the rescaling, the results (solid lines in Fig. [S] (a) 
and (b)) exhibits better agreement with the DFT results. 
Thus despite the simplicity of the TB model, it describes 
the various aspects of the behavior of the polarization 
well including the width dependence. 

To experimentally measure this proposed effect, one 
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FIG. 5. (Color) Comparison between the results of DFT cal¬ 
culations and those with the TB model, (a), (b) Polarization 
obtained by DFT calculations (squares) compared with re¬ 
sults with the TB model (lines). The results of the TB with 
(without) rescaling by the gap size are shown as a solid line 
(dotted line) (see text), (a) {N = 9,12,15) are for N = 31, 
and (b) {N = 7,10,13) are for N = 3Z-I-1. Note that it is for a 
spinless system, and the results here should be multiplied by 
two to compare with experiments, (c) Induced polarization 
P by DFT calculation divided by {N -|- 1)^. (d) Relationship 
between out-of-plane electric field |E| versus on-site potential 
energy difference U. U is scaled by the hopping amplitude 
t = 2.6eV 0]. 


needs a bilayer nanoribbon with well-defined edges and 
width. For single-layer graphene nanoribbons, well- 
defined edge orientations have been demonstrated [IB- 
[ 2 ^, and it might be realized also for bilayer graphene. 
For the bilayer graphene, interlayer electric field up to 
0.3V/A has been achieved [^. and therefore the pro¬ 
posed effect with polarization up to ^ —0.12e per spin 
for N = 15 is expected to be realizable experimentally. 
We also have calculated the effect of periodic modula¬ 
tions of the width to check the edge disorder effect via 
supercell approach and confirmed that t he p olarization 
survives the weak modulations considered ll| . Neverthe¬ 
less, since the effect is sensitive to the ribbon width, the 
proposed effect will disappear in the presence of strong 
disorder. We note here that the in-plane polarization by 
an interlayer bias can be expected for a wide variety of 
atomic-layer compounds, as long as the symmetry crite¬ 
rion for its emergence is satished. As an example, a bi¬ 
layer armchair ribbon of transition metal dichalcogenides 
in the 2H stacking satishes this criteria. Moreover, our 
calculation show induced polarization in AA’-stacked bi¬ 
layer boron nitride nanoribbons 11| . Such a wide choice 
of candidate materials provides us with many chances for 
experimental verifications of our theory. 

To conclude, we theoretically show that the AB- 


stacked graphene nanoribbon with armchair edges has a 
polarization along the ribbon direction, when interlayer 
bias voltage is applied. This is shown both by the simple 
tight-binding model and the ab initio calculations. In 
particular, the linear response to the interlayer voltage 
shows different signs for the cases N = 31 and N = 3Z-I-1, 
and it is fully understood by means of a simple two-band 
model. 
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MEXT, Japan (No. 26287062, No. 25107005 and No. 
25104711), JSPS Research Fellowships for Young Scien¬ 
tists, and MEXT Elements Strategy Initiative to Form 
Core Research Center (TIES). 
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FORMULA OF THE POLARIZATION IN TERMS OF THE BLOCH WAVEFUNCTIONS 


In our calculation of the polarization P in terms of the Bloch wavefunctions, we used the formula 

OCC. Q 

dk'^{ukn \ \ukn) , 

n 

where \ukn) is the Bloch wavefunction satisfying the cell-periodic gauge condition 

Uk,n{r) = 



(SI) 


(S2) 


and the summation is taken over the occupied states below the Fermi energy. Here k is the Bloch wavenumber and 
G = 27r/a is a reciprocal lattice vector. In numerical calculation, the differentiation in terms of k in Eq. (ISII) should 
be replaced by a difference in k. Such a formula with this replacement is discussed in detail in Ref. 0, and we followed 
this formalism for the calculation of the polarization. 


POLARIZATION CALCULATED FROM THE TIGHT-BINDING MODEL FOR VARIOUS WIDTHS 


We show numerical results of the polarization induced by the interlayer bias for graphene nanoribbons (GNRs) with 
various widths iV, calculated from the tight-binding model. Some results are shown in Fig. 2 in the main text, and 
we show more examples for various N in Fig. ISll For N = 31 the results are shown in Fig. ISlh . c {N = 9, 15, 21), and 
in Fig. 2a, b {N = 27, 39). For V = 3/ -I- 1 the results are shown in Fig. ISlb . d {N = 10, 16, 22), and in Fig. 2c, d 
{N = 28, 40). In the wide range of U iFi.g lSTI) . the polarization P oscillates as a function of U, which is attributed to 
crossings of minigaps. There are more oscillations for larger N, which reflects the fact that there are a larger number 
of minibands for wider ribbons. On the other hand, the oscillation amplitude becomes gradually smaller for wider 
ribbons (see Fig. 2) because the contribution from each miniband becomes relatively smaller. On the other hand, in 
the regime U ^ 0, we showed from the two-band model that the polarization is linear in U with its slope scales with 
(iV -I- 1)^. This is roughly reproduced for the results shown in Fig. ISlb and d. 




FIG. SI. Polarization induced by the interlayer bias U for GNRs with various widths N, calculated from the 
tight-binding model. It is shown as a function of the interlayer bias U. a and c are for the class N = 31 and b and d are 
for the class N = 31 + 1. 

It may look strange that for a large N limit, i.e. the 2D graphene limit, the polarization has a different asymptotics 
for N = 31 and N = 31 + 1. It is in fact reasonable; the polarization per area is proportional to P divided by the 
width, and therefore goes to zero for V —> oo. 
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EIGENVALUES AND EIGENVECTORS OF THE TIGHT-BINDING MODEL 


Here we derive the two-band low-energy effective Hamiltonian from the eigenstates of the tight-binding model of 
the bilayer GNR with armchair edge, without the interlayer bias voltage: 


Ht — tijc\cj. (S3) 

<i,0> 

Firstly, to obtain the two-band Hamiltonian, we diagonalize the tight-binding model at fc = 0 according to Ref. [l|. We 
set the eigenvector of the Hamiltonian 

N N N N 

= 0)) =Y^ am \mA) +^hm\mB) + ^ \mA') + ^ 1™^'). (S4) 

m —1 m —1 m —1 m —1 

where A and B represent sublattices in the lower layer, and A' and B' in the upper layer. Here, Om, &m, a'm, and b'^ 
are expansion coefficients. From the tight-binding model at fc = 0, we obtain 


— ^(^m—1 “b “b ^m+l\ 

~ “b CLrn “b Clm+l) “b (S6) 

+t±bm, (S7) 

>^b'm = t{a'm-i+a'm + am+i), (S8) 


where e represents the energy eigenvalue. To solve the above equations (IS5ll - (IS8ll . we introduce new coefficients 


Then, we obtain 

£ra± = t(/3^_i-b/3^-b/3^+i), (SIO) 

£/3m = i(am-i+a™ + a^+i)±t_L/3^. (Sll) 

Since = 0 = (3^, the solutions have the form cx sin(m0) and (3^ B^sm.(m0), where 0 is a constant 

(0 < 0 < tt). We rewrite the equations (ISlOll and (ISllI) in the matrix form 


f 0 t(2cos6l-Hl)\ _ ±fA^\ 

l^t(2cos6»-bl) ±t_L )\B^)~^ 
and its eigenvalues are obtained analytically. 

= ±Y + {y) +^^(cos 6'-b 1)2, 


(S12) 


(S13) 


where q = ±1. From the boundary condition, the coefficients must vanish when m = 0 and V -|- 1 and we get 


dr = ^ = 1>2,...,V. (S14) 

Thus, we get AN eigenvalues and eigenvectors, 

±Y Tldr, 

1 f T^Y+ddr \ 

^2dr{dr T g^) v(2cos6»r-b l)y 

We put dr = \J (t_L/2)^ -I- t^(2 cos^r + 1)^ for notational simplicity. These energy eigenvalues can become zero only 
when cosdr = —1/2. When N = 3/ -I- 2 it can be satisfied for r = 21 + 2, and the energy bands are gapless at 
fc = 0 because e2i’3^2 — ^ 21+2 — 1^® other hand, when N = M or 31 + 1, cos 9r = —1/2 cannot be satisfied the 

energy bands are gapped. Therefore, the polarization Py(U) can be defined in the nanoribbons with widths N = 31 
or 3Z -I- 1. In these cases, the energy eigenvalues closest to zero are and and the corresponding eigenstates 


(515) 

(516) 











7 


1+) = |r = 2Z + 1, —, +) and |—) = |r = 2/ + l,+,— ) considerably contribute to the polarization. For brevity, we write 
±50 = and 


50 =—^ + d2i+i, (S17) 

l±) = / = X 

V2Em=iSin^(™^'2;+i) 

N 

X] sin(TO6»2/+i)(A|:fi(|myl) =F \mB')) + T 1™^')))- (S18) 

m—1 

Hereafter, we omit subscripts 2/ + 1 except for that in O 21 + 1 ] for example, we write d = d 2 i+i- 


CONSTRUCTION OF THE TWO-BAND HAMILOTNIAN 


To elucidate the polarization in the weak interlayer bias voltage, we construct a two-band low-energy effective 
Hamiltonian for the space spanned by |±). Therefore, we add the interlayer bias voltage U to the tight-binding model 
as a perturbation; 




2 cj Cj , 


(S19) 


We retain terms up to the linear order in k and the nonzero matrix elements to this order are given by 
{mA\ Ht \nB) ={mA'\ Ht \nB') 

=t{S 

m,n —1 -k s m,n ± ) ± 2 (^m,n 2 *^^+1,21)5 


{mA'\Ht \nB) =t±Sjn,n, 




(mA\ Hu \nA) ={mB\Hu \nB) = -—5m,n, 
{mA'\Hu \nA') =(mH'| Hu \nB') = 
Therefore, the nonzero matrix elements of Ht and Hu are 

(±1 Ht |±) =±5o. 


(-|-| Ht |-) =ikf — {l - cos52i-i-i)(2cos52Z+i ± 1), 


3d 


{+\Hu\-) = 


Ut^ 


Hence, we obtain the two-band Hamiltonian HeS = Ht + Hu, 

Heff{k,U) = hiUax + h2kcFy + goa^, 

where 

hi = h 2 = - cos52Z+i)(2cos52Z-hi ± !)■ 

The eigenvalues and eigenstates of this two-band Hamiltonian are 


Ek,±= ±\j[hiUy + (^2^)2 -H = ±g, 


|wfe,±) = 


50 ±5 


x/ 25 ( 5 ± 5 o) \hiU + ih2kJ ■ 


(520) 

(521) 

(522) 

(523) 


(524) 

(525) 

(526) 


(S27) 


(S28) 


(529) 

(530) 








POLARIZATION FROM THE TWO-BAND HAMILTONIAN 


We focus on the region of the weak interlayer bias voltage to clarify the difference of the slope of the polarization 
at {7 ^ 0 for two classes N = 31 and TV = 3^ -I- 1. The polarization P is given by 


p= r du'j{u'), 

Jo 


where 


occ. unocc. 


/“ {Ukn\ \Ukm) {Ukm\ gif \Ukn) , 

Z. {Ekn-Ekm? 


(S31) 


(S32) 


Here, j{U) = dP/dU, and n and m are band indices for occupied bands and unoccupied bands, respectively. Ekj 
and \ukj) are the jth eigenvalue and the eigenstate of the Hamiltonian. In the present case n = — and m = -I-, and 
therefore P and j{U) are calculated by using the two-band Hamiltonian derived in section „ 


- , e / hih 2 UTr/a 

nU) = -^^arctan - , , . 


27r \goV9o + {hiUy + (/i27r/a)2 


3{U) = --h^h2- 


50 


2 a {{hiUY -f gl) ^{hiUy +gl + 
In particular, the slope of the polarization at t/ = 0, i.e. j(0), is given by 


j( 0 ) = -7r^i^2- 


1 


2a ffoVso + (/i27r/a)2 
Because t±_ > 0, hi and go are positive, and we have 

'2 


Here, 


sgn(j( 0 )) = -sgn(Ti 2 ) = sgn(^- 7 r - 6 » 2 /+i) . 


In- 921+1 = 1 

3 1 N = 3l + 1 


. 3(Af-l-l) 


Therefore, the sign of j(0) = i9P/cI?7|£/=o is given by 


dP 


I < 0 N = 3l 
c/=o I > 0 = 3/ -b 1 


(533) 

(534) 

(535) 

(536) 

(537) 

(538) 


which well agrees with the numerical results of the tight-binding model. Asymptotic behavior of j(0) (Eq. (IS35I) ') is 
evaluated for large TV, where go <C /i 27 r/a; 


j( 0 ) = -sgn(/i 2 ) 


t Le 


^le ^ ^_ 

27rgo 8TTd{d— 


where ± is — for TV = 31 and -b for TV = 3Z -b 1. By using equation (IS37I) it is approximated as 


(S39) 


(S40) 


EFFECT OF PERIODIC MODULATION OF THE RIBBON WIDTH 


We numerically calculate the polarization in bilayer GNRs with weak periodic modulations of the ribbon widths by 
using the tight-binding model when the interlayer bias voltage is weak. We consider two cases of periodic modulations 
by changing the width of each layer in various ways. 
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Firstly, we discuss an effect of the difference of the widths of the upper and lower layers. Figure [S2] shows the 
polarization in the armchair bilayer GNRs when the upper and lower layers have the different widths. Then, we 
calculate the polarization by changing the width for the upper layer Nu while fixing that for the lower layer TVl, 
except that TVl = 2 or iVy = 2 (mod 3) since the energy bands are gapless at U = 0. When the system is 

the perfectly stacked bilayer armchair GNRs in Fig. 1. The results in Fig.[S2]b and c correspond to the polarization 
for = 12 = 0 and A^l = 13 = 1 (mod 3), respectively. Gonsequently, we can see that all the slopes of polarization 
in Fig. [S2] have the same sign regardless of the width of the upper layer Ntj. Furthermore, even if the bilayer is 
composed of two layers with the width 0 and 1 (mod 3), we find that the sign of the slope is equal to that of the 
perfectly stacked bilayer GNRs with the narrower width although the magnitude becomes small. Therefore, when the 
widths of the upper and lower layers are different, the polarization behaves like the perfectly stacked armchair GNRs 
with the width min(N^Ui Nl). 



FIG. S2. Polarization in the bilayer composed of two layers with different widths, a shows structure of the GNRs 
with the different upper and lower layer. The unit cell is described by the red dashed-line box. b and c show the polarization 
for Nu = 12 and Nu = 13, respectively, b is for Nu = 31, and c is for Nu = 31 -f 1. 


Secondly, we calculate the polarization of bilayer GNRs when the widths of the upper and lower layers alternates 
between two values Ni and N 2 {Ni > N 2 ), as shown in Fig. [S3] a having no dangling bonds. In this system, the 
primitive translation vector doubles. In particular, when Ni = N 2 , the system corresponds to the armchair bilayer 
GNRs. In this case, when we calculate the polarization, we change only Ni and fix N 2 . The results of the polarization 
for various widths {Ni,N 2 ) are shown in Fig. IS3I Figure |S3| b and c show the polarization for N 2 = 12 = 0 and 
N 2 = 13 = 1 (mod 3), respectively. As a result, as Ni becomes larger, we find that the magnitudes of the polarization 
for N 2 = 0 (mod 3) are enhanced, while those for N 2 = 0 (mod 3) are suppressed. Nevertheless, the sign of the slope 
of the polarization is unchanged from that of the perfect armchair GNRs with the width Ni = N 2 . Therefore, we 
can see that the polarization in the weak interlayer bias voltage is dominated by the narrow part of GNRs, which is 
similar to the previous case. 

From the above results, we find that small variations of the width do not affect the sign of the slope of the polarization 
in the weak interlayer bias voltage. In other words, even though the edges of the bilayer graphene nanoribbons are 
not completely perfect, the nontrivial dependence of the polarization on the width appear like the perfectly stacked 
bilayer GNRs with armchair edges. Thus, the in-plane polarization in response to the interlayer voltage survives even 
if the edges have weak periodic modulation of the ribbon width, as long as the energy bands are gapped. Nevertheless, 
since the effect is sensitive to ribbon width, the proposed effect will disappear in the presence of strong disorder. 
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u/t 


u/t 


FIG. S3. Polarization in the bilayer consisting of the stacked monolayer with two widths, a shows the structure 
of the bilayer nanoribbons with two widths. The red dashed-line box describes the unit cell, b and c are examples of the 
polarization for A^u = 12 and = 13, respectively, b is for A'^u = 31, and c is for A% = 3Z -I- 1. 


POLARIZATION IN THE BILAYER BN NANORIBBON 


We explained the emergence of the polarization in bilayer graphene nanoribbons by symmetry argument. Therefore, 
other nanoribbons of atomic-layer compounds can have a finite polarization along the ribbon direction in response 
to the interlayer voltage, when the symmetry criterion for such a response is satisfied. To confirm this, we compute 
the polarization of hydrogen terminated bilayer BN nanoribbons from hrst-principles calculations. Here, we consider 
so-called AA’-stacked bilayer BN nanoribbons with the armchair edges and the geometry is fully optimized (see Fig.lSdl 

(a) ). As in the case of the bilayer armchair GNRs, xz-mirror symmetries are broken in this structure. When the 
interlayer bias voltage is zero, inversion symmetry is preserved and the polarization is zero. The interlayer voltage 
breaks the inversion symmetry, leading to nonzero polarization along the ribbon, as we see in the following. Figure [S4] 

(b) shows calculated in-plane polarization as a function of the interlayer voltage. As we can see, a finite polarization 
appears as expected. Note that the size of the polarization is rather small compared to those in the GNRs. One 
reason is that the band gaps of these nanoribbons are relatively large 4.3 eV both for A = 7 and N = 9 at U = 0). 
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Electric field (V/A) 


FIG. S4. Polarization in the BN nanoribbons, a Top view of the hydrogen terminated bilayer BN nanoribbon with the 
armchair edges for N = 7. B (N) atoms on the lower layer are just below the N (B) atoms on the upper layer, b Electric field 
dependence of the polarization obtained by DFT calculations for bilayer BN nanoribbons with the armchair edges. 




